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Abstract 

Background: In East Asia, an increasing number of studies on temperate forest tree species find evidence for 
migration and gene excliange across tlie East Cliina Sea (ECS) land bridge up until the last glacial maximum (LGM). 
However, it is less clear when and how lineages diverged in this region, whether in full isolation or in the face of 
post-divergence gene flow. Here, we investigate the effects of Quaternary changes in climate and sea level on the 
evolutionary and demographic history of Platycrater arguta, a rare temperate understorey shrub with disjunct 
distributions in East China (var. sinensis) and South Japan (var. arguta). Molecular data were obtained from 14 
P. arguta populations to infer current patterns of molecular structure and diversity in relation to past (Last 
Interglacial and Last Glacial Maximum) and present distributions based on ecological niche modelling (ENM). 
A coalescent-based isolation-with-migration (IM) model was used to estimate lineage divergence times and 
population demographic parameters. 

Results: Combining information from nuclear/chloroplast sequence data with nuclear microsatellites, our IM 
analyses identify the two varieties as genetically distinct units that evolved in strict allopatry since the mid- 
Pleistocene, c 0.89 (0.51-1.2) Ma. Together with Bayesian Skyeline Plots, our data further suggest that both lineages 
experienced post-divergence demographic growth, followed by refugial isolation, divergence, and in the case of 
var. arguta post-glacial admixture. However, past species distribution modelling indicates that the species' overall 
distribution has not greatly changed over the last glacial cycles. 

Conclusions: Our findings highlight the important influence of ancient sea-level changes on the diversification of 
East Asia's temperate flora. Implicitly, they challenge the notion of general temperate forest expansion across the 
ECS land bridge, demonstrating instead its 'filter' effect owing to an unsuitable environment for certain species and 
their biological (e.g., recruitment) properties. 
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Background 

Phytogeographers have long recognised the Sino-Japanese 
Floristic Region (SJFR) of East Asia as the worlds major 
centre of temperate plant diversity and endemism [1]. 
Much of the potential primary vegetation of this vast re- 
gion is composed of warm-temperate deciduous (WTD) 
forest, as presently found scattered in mid-elevation 
subtropical (Central/South/East) China, predominant in 
low-elevation North China and the Korean Peninsula, 
and disjunctively distributed in the main islands of Japan 
[1,2]. Fossil pollen analyses have previously indicated that 
during the Last Glacial Maximum [LGM: c. 21,000- 
18,000 yr before present (BP)], the habitat of East Asian 
WTD forests in the northern parts of their range (e.g., in 
North China and North Japan) contracted, mainly in re- 
sponse to increased aridification [3,4]. However, palaeo- 
biome reconstructions suggest that these forests also 
expanded across the large expanses of continental shelf 
(c. 1 million km^) that emerged in the East China Sea 
(ECS) as a consequence of sea level lowering by c. 85- 
130/140 m during cold periods [5-7]. In consequence, it 
is now widely believed that a near continuous belt of 
WTD forests spanned the ECS continental shelf during 
the LGM (and possibly earlier cold periods), connecting 
presently disjunct populations in East China, South 
Japan, and the Korean Peninsula [2,8]. If correct, this 
land bridge may have served as a 'dispersal corridor' 
[9,10], allowing intermittent migration of most WTD 
forest-restricted plant species from the Asian mainland 
into Japan (or vice versa), and/or periodic secondary con- 
tact and gene flow among formerly isolated populations, 
possibly up until the last shelf submergence (c. 16,000- 
7,000 yrBP; [11]). However, there is accumulating but still 
scanty evidence that the ECS land bridge also acted as a 
'filter' in selectively hampering or even preventing the 
dispersal of certain plant species, while allowing those 
able to tolerate the environmental conditions of this 
palaeo-landscape to disperse more freely (see below). 

Support for the 'recent dispersal corridor hypothesis' 
comes from phylogeographic data of two wide-ranging 
and abundant WTD forest tree species of the SJFR 
{Cercidiphyllum japonicum [12] and Kalopanax septem- 
lobus [13]). In both instances, there is a near lack of differ- 
entiation for chloroplast and/or nuclear DNA sequences 
across the ECS, consistent with ecological niche modelling 
(ENM) predicting large expanses of suitable habitat for 
each species on the ECS land bridge during the LGM. By 
contrast, a high level of genetic differentiation across the 
ECS has been identified in several rare understorey herbs 
and shrubs restricted to the mountainous WTD forests 
of East China and South Japan {Croomia japonica [14], 
Platy crater argutal Kirengeshoma palmata [15,16] and 
Ligularia hodgsonii [17]). In all of these latter cases, the 
estimated times of trans-ECS divergence based on 



chloroplast (cp) DNA fall into the early-to-mid Pleisto- 
cene, suggesting that the ECS land bridge imposed a 
formidable barrier to dispersal during the last glacial 
cycle(s) (reviewed in [18]). However, all of these previous 
time estimates relied on a single locus, i.e., chloroplast 
(cp) DNA, which thus renders them subject to error from 
the inherent stochastic nature of the coalescent [19-21]. 

The aim of the present paper is to re-examine the 
evolutionary history of one of these understorey shrubs, 
Platycrater arguta Siebold & Zucc. (Hydrangeaceae), using 
a multi-locus approach. Platycrater arguta, the only 
species of this monotypic genus, is a small deciduous 
shrub endemic to the mountainous WTD forests of East 
China and South Japan, where respectively two varieties, 
var. sinensis and var. arguta, are recognized [22-24]. In 
the cpDNA study previously conducted [15], these taxa 
were found to comprise distinct phylogroups, whose 
likely vicariant origin was dated to the mid-Pleistocene 
(c. 0.89 Ma). Here we use a broader sampling of P. arguta, 
and present additional datasets of two nuclear DNA 
(nDNA) sequence markers (ITS, Tpi) and nuclear micro- 
satellite (nSSR) loci to infer a more robust divergence and 
demographic history of this species. Specifically, our goals 
were (i) to provide a refined time-scale for the divergence 
of var. sinensis and var. arguta by fitting all four datasets 
(ITS, Tpi, nSSRs, cpDNA) to an 'isolation with migration' 
(IM) model [20,25]; (ii) to quantify the amount of post- 
divergence gene flow between them while accounting for 
potential changes in effective population sizes; and (iii) to 
model the species' potential distributions in response to 
past climatic changes, specifically from the Last Intergla- 
cial (LIG/Eemian: c. 130,000-114,000 yrBP) [26] through 
the LGM to the present day. We attempted to reconcile 
these distribution patterns with our genetic data to deter- 
mine the role of the exposed ECS shelf as a 'corridor' or 
'filter' during lineage divergence at the time of the last 
glacial cycle (s). Finally, together with ENM, the added 
value of three nuclear data sets allowed for more confi- 
dence in the interpretation of the supposedly contrasting 
population histories of each variety than was derived pre- 
viously from a single dataset (cpDNA) [15]. 

Methods 

Plant material and sampling design 

We sampled seven populations of each of var. sinensis 
and var. arguta, {212 individuals in total) covering the 
species' entire range in East China and South Japan 
(Additional file 1: Table SI). Ten of the populations (100 
individuals) were analysed previously for cpDNA [15], 
whUe four were newly sampled (C3-C5 and J2). As 
P, arguta is nested within a paraphyletic Hydrangea 
L., but with yet undefined sister species [27,28], we ar- 
bitrarily designated H, chinensis Maxim., collected from 
Yandang Mt. (China), as one of the outgroups together with 



Qi et al. BMC Evolutionary Biology 2014, 14:41 
http://www.bionnedcentral.conn/1471-2148/14/41 



Page 3 of 16 



other members of Hydrangeaceae (see below). Voucher 
specimens of this species and all sampled populations 
of P. arguta are stored at the Herbarium of Zhejiang 
University (HZU). 

Molecular protocols 

Total genomic DNA was extracted from silica-dried leaf 
tissue by the cetyltrimethyl ammonium bromide (CTAB) 
method [29]. The entire internal transcribed spacer (ITS) 
region of nuclear ribosomal (nr) DNA was sequenced in 
72 individuals of P, arguta (4 to 8 individuals per popula- 
tion; Additional file 1: Table SI) and one individual of 
H, chinensis, following the methods described in [14]. 

After preliminary screening of six low-copy number nu- 
clear genes using primers developed by Strand et al. [30] 
(ADHX2F-4R, CAMX1F-2R, CHIX1F-4R, CHSXIF- 
2RN, GPDX7F-9R, TPIX4FN-6RN), we chose the triose- 
phosphate isomerase {Tpi) gene (TPIX4FN-6RN) for the 
full survey of 173 individuals (5-20 individuals per 
population, plus H, chinensis; Additional file 1: Table SI) 
because this gene region was single copy in most cases 
and proved to be sufficiently polymorphic [31]. Direct 
sequencing of Tpi was carried out using the same pro- 
cedure as described in Zou et al, [32]. Chromatograms 
of ITS and Tpi with additive peaks were further ana- 
lysed by inferring the identity of the two haplotypes 
within a heterozygote through haplotype subtraction 
[33,34]. When the chromatogram quality did not permit 
this procedure, PGR products were cloned using a 
pMD18-T vector system (Takara Biotechnology, Dalian, 
Liaoning, China) according to the manufacturer s proto- 
col. Six to ten clones were sequenced per individual 
when sequences contained two or more polymorphic sites. 
The two sequences of a heterozygote were separated by 
comparing sequences of the PGR product and the cloned 
sequences. All haplotype sequences have been deposited 
in GenBank [accession numbers for P, arguta: JQ978221- 
JQ978253 (ITS) and JQ978254-JQ978284 {Tpi); for H, 
chinensis: KF559183 (ITS) and KG853063 {Tpi)], 

All 272 P, arguta individuals were genotyped at seven 
nuclear dinucleotide-repeat microsatellite loci (Pal-Pa3, 
Pa5-Pa8) (Additional file 2: Table S2) according to the 
methods described in Qi et al, [35]. Fluorescently labelled 
PGR products (HEX or 6-FAM; Applied Biosystems, Foster 
Gity, GA, USA) were separated on a MegaBAGE 1000 
(GE Healthcare Biosciences, Pittsburgh, PA). Alleles were 
scored manually with the aid of GENETIG PROFILER 
2.2 (GE Healthcare Biosciences) using the ET550-R size 
standard. 

Nuclear DNA sequence analyses 

Sequences of ITS and Tpi were aligned and edited in 
GENEIOUS 4.8.2 [36]. For each gene region, we esti- 
mated haplotype diversity {h) and nucleotide diversity 



(tt) for each population, for each variety, and for the 
whole data set. Tests for non-neutral evolution were 
conducted by computing Tajimas D [37] and Fu & Lis 
[38]. In addition, we estimated Fus Fs [39] and 
Ramos-Onsins and Rozas' R2 [40] to detect population 
growth. Gritical values for these statistics were obtained 
using 10^ coalescent simulations. Recombination was 
calculated as the minimum number of recombination 
events {Rm) using the four-gamete test [41]. All of the 
above analyses were performed in DNASP 4.1 [42] and 
ARLEQUIN 3.11 [43]. 

Phylogenetic relationships of ITS and Tpi haplotype 
sequences were inferred using maximum parsimony (MP) 
and maximum likelihood (ML), with gaps (indels) treated 
as missing data. Heuristic MP searches were performed in 
PAUP* 4.0bl0 [44] using the same settings as described in 
Qiu et al, [16]. The ML analysis was conducted under the 
GTR + r substitution model using RAXML 7.2.8 [45]. 
Node support was assessed using 1,000 'fast bootstrap' 
replicates. The following species were chosen as out- 
groups: Hydrangea chinensis (ITS, Tpi) as well as H, 
anomala D. Don and Schizophragma hydrangeoides 
Siebold & Zucc. (only ITS; downloaded as GenBank 
accessions JF976651 and GU98303, respectively). The 
alignments and phylogenetic trees were deposited in 
TreeBASE (submission number 15354). In addition, we 
constructed haplotype networks for each dataset using 
the 95% statistical parsimony criterion implemented in 
TGS 1.21 [46]. For ITS, however, we had to increase 
the TGS connection limit to 50 steps to link the diver- 
gent networks of the two varieties. Indels were treated 
as single mutation events, and coded as substitutions 
(A or T). Population differentiation for unordered (Gst) 
and ordered (A/st) haplotypes were obtained with PERMUT 
to test whether A^st is significantly larger than Gst 
(1,000 random permutations), indicating the presence 
of phylogeographic structure [47]. Analyses of molecu- 
lar variance (AMOVAs) were carried out in ARLEQUIN 
using F-statistics. Sequence variation was hierarchically 
partitioned between the two varieties, among popula- 
tions within varieties, and within populations. The sig- 
nificance of all estimated fixation indices was tested 
using 10,000 permutations as described in Excoffier 
etal, [48]. 

Nuclear microsatellite analyses 

Measures of nSSR diversity were assessed for each popu- 
lation, and across all loci, by calculating in FSTAT 2.9.3 
[49] the total number of alleles (A/a), allelic richness {R^) 
standardized for 5 individuals using rarefaction), expected 
gene diversity (//s), and the inbreeding coefficient (Fjs). 
Differentiation between populations was computed in 
ARLEQUIN using i?sT [50], which assumes a stepwise 
mutation model (SMM; [51]), and thus may be more 
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realistic for microsatellite data than other measures (e.g., 
Fst) based on the infinite alleles model (I AM). 

Overall population structure was examined using the 
Bayesian clustering algorithm implemented in STRUC- 
TURE 2.3 [52]. This program was run 10 times on individ- 
ual multi-locus nSSR genotypes for a number of clusters 
(/<), ranging from 1 to 14 (the number of localities), using 
a burn-in length of 50,000 and run length of 500,000 itera- 
tions. We used the admixture model without prior infor- 
mation on sample population membership and allowed 
allele frequencies to be correlated among clusters [53]. 
We plotted the posterior probability of the data [lrLP(£))] 
and the ad hoc statistic AK [54] for determining the most 
likely K 

Hierarchical AMOVA was performed using 7?-statistics, 
and partitioned as described above. In order to detect 
genetic signatures of recent population bottlenecks in 
the nSSR dataset, we applied two tests implemented in 
BOTTLENECK 1.2.02 [55]: (i) Wilcoxons sign-rank test 
for heterozygote excess; and (ii) the 'mode-shift test' for 
detecting a shifted, rather than an equilibrium, L-shaped 
distribution of alleles. 

Divergence time, effective population sizes, and gene 
flow 

We used the 'isolation with migration' (IM) coalescent 
model as implemented in IMA [25] to estimate population 
rate parameters (0) and effective population sizes (A/e) of 
var. sinensis in China (0c)> var. arguta in Japan (0j) and 
their common ancestral population (©a; 0 = 4A/eu), as 
well as bidirectional migration rates (mc-j and mj_c; M = 
m/u) and divergence times (r = ^u) between the two var- 
ieties. All parameters in the IM model are scaled by the 
neutral mutation rate (u). For this analysis we jointly 
employed the present nuclear datasets (ITS, Tpiy nSSRs) 
and previously obtained sequences of two cpDNA regions 
{psbA-trnH, trnD-trnE; [15]). Since IMA assumes no 
recombination within loci, the nDNA sequence data 
were analysed by using only maximally informative and 
nonrecombining blocks per individual as identified by 
the program IMGC [56]. 

Both nDNA and cpDNA sequence data were analysed 
under the HKY model rather than the infinite sites model 
(ISM) as we regularly found some sites with more than 
two changes in these data sets (note that HKY and ISM 
are the only models of nucleotide evolution implemented 
in IMA, whereby only HKY allows multiple changes 
at a site). Corresponding values of u were calculated 
as u = //kg, where ^ is the number of substitutions 
per site per year, k is the average sequence length 
under study in base pairs (excluding indels), and g is 
the generation time in years (i.e., age of first reproduction), 
which is five years in P. arguta as observed under cultiva- 
tion [15]. As substitution rates for this species are 



unknown, we assumed the following mean values and 
confidence intervals (CI) in substitutions per site per year: 
ITS (for woody perennials), 2.15 x 10"^ (0.38-7.83 x 10"^) 
[57]; Tpi, 7x10"^ (0.5-3.0x10"^) [58,59]; and cpDNA, 
1.52x10"^ (1.0-3.0x10"^) [58,60]. The nSSR loci were 
assumed to fit a stepwise model of mutation (SMM), with 
a mean mutation rate ((i) of 10"^ mutations per locus per 
year. Although the assumptions underlying this mutation 
rate are debatable (e.g., [61,62]), this rate falls close to the 
median of average values (c. 7.7 x 10"^) reported for a 
wider range of plant species [63-67], and has also been re- 
cently employed for other woody perennials [68,69]. The 
geometric average mutation rate of the four marker sets 
was used to rescale the IMA parameter estimates from the 
combined analysis. The inheritance scalars were set to 1 
for the nuclear markers, and 0.5 for cpDNA, the latter 
value reflecting the expected effective population size (A/g) 
of maternally inherited cpDNA in a hermaphroditic 
plant relative to an autosomal locus [70] (note that ma- 
ternal inheritance of cpDNA has been demonstrated in 
Hydrangeaceae [71]). Markov chain Monte Carlo (MCMC) 
simulations were conducted for 10^ steps by using a linear 
heating scheme and 10 Metropolis-coupled chains with a 
burn-in period of 10^ steps. To verify convergence on the 
same parameter values, we ran this analysis three times 
with different random seeds. Only estimates whose poster- 
ior distribution dropped to zero within the prior intervals 
were trusted. 



Historical demography 

To further examine the historical demography of each 
variety, we estimated the shape of the population growth 
function through time by constructing Extended Bayesian 
Slcyline Plots (EBSPs) as implemented in BEAST 1.7.5 [72]. 
This coalescent-based, nonparametric Bayesian MCMC 
algorithm incorporates multi-locus data to reduce esti- 
mate errors associated with single genes and increases 
the power to detect demographic dynamics [73]. Analyses 
were performed for each variety across the sequence data 
sets (ITS, Tpiy cpDNA) assuming a strict molecular clock. 
The mutation rate priors for each locus were identical 
to those used for IMA. Based on the Akaike Informa- 
tion Criterion (AIC) as implemented in JMODELTEST 
[74], we selected the GTR + F and HKY substitution 
models for the nuclear and cpDNA sequences, respect- 
ively. Bayesian MCMC chains were run for 10 million 
generations, with a sampling frequency of every 1,000 
generations, whereby the first 5,000 samples were dis- 
carded as burn-in. Convergence of the MCMC chains 
was inspected using TRACER 1.5 [75] by visually 
checking that effective sample size (ESS) for all rele- 
vant parameters were well above 200. Skyline plots 
were visualized using EXCEL. 
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Present and past distribution modelling 

To infer distributional changes of P. arguta since the last 
interglacial, we produced ENMs for three time periods 
(present, LGM, LIG) using the maximum entropy method 
implemented in MAXENT 3.2.1 [76]. In addition to the 
distribution records included in this study, records were 
sourced from the Chinese Virtual Herbarium (www.cvh. 
org.cn), the National Specimen Information Infrastructure 
of China (www.nsii.org.cn), the Kyoto University Herbar- 
ium (KYO), and Red Data Books for the Aichi and Mie 
Prefectures of Japan [77,78]. Based on a total of 59 records 
(China: 20; Japan: 39), a current distribution model was 
developed using six bioclimatic data layers (biol, 12, 16, 
17, 18, 19; WorldClim dataset, [79,13]) at 2.5 arc-min 
resolution. This model was then projected onto the set of 
climatic variables simulated by the Model for Interdiscip- 
linary Research on Climate (MIROC) 3.2 [80] to infer the 
extent of suitable habitat during the LGM and the LIG 
[81-83,13]. The accuracy of each model prediction was 
then tested by calculating the area under the 'Receiver 
Operating Characteristic (ROC) Curve' (AUC; [84]). We 
note that the ENM for the LGM explicitly accounted for 
changes in palaeo-coastlines (-110 m than at present) and 
palaeo-climate surfaces on the exposed ECS continental 
shelf [7,85]. 

Results 

Nuclear sequence characteristics 

The ITS sequences of the 72 P, arguta individuals (14 
populations) from East China (34) and South Japan (38) 
were aligned with a total length of 737 base pairs (bp), 
revealing 36 nucleotide substitutions and three 1-bp indels. 
Together, these 39 polymorphic sites yielded 33 ITS hap- 
lotypes (ribotypes; H1-H33) (Additional file 3: Table S3). 
For the Tpi locus, only one or two distinct sequences were 
detected in each of the 173 individuals surveyed, 28 of 
which were found to be heterozygotes. In total, 31 haplo- 
types (T1-T31), ranging from 309 bp to 314 bp, were des- 
ignated based on 39 substitutions and three small (< 5 bp) 
indels (Additional file 4: Table S4). None of the loci 
showed significant deviation from neutral expectations for 
Tajimas D or Fu and Lis at the species or variety levels 
(Additional file 5: Table S5). Demographic tests of Fus Fs 
revealed negative but nonsignificant values, while R2 was 
consistently positive and significant, suggesting demo- 
graphic growth (Additional file 5: Table S5). The mini- 
mum number of recombination events {Rm) estimated for 
ITS and Tpi were 15 and 2, respectively. 

Genetic diversity, haplotype relationships and genetic 
structure at ITS and Tpi 

There were high and comparable levels of total haplotype 
(/zt) and nucleotide (ttt) diversity at the species-level (ITS/ 
Tpi: hr = 0.96/0.92; ttt = 0.03923/0.02133), and the same 



was found for each variety (region), whereby each of the 
seven populations from East China {ITS/ Tpi: = 0.88/ 
0.90; ttt = 0.0151/0.0108) and South Japan {ITS/Tpi: = 
0.95/0.87; ttt = 0.0188/0.0139) harboured broadly similar 
levels of diversity (Additional file 1: Table SI). Concordant 
with the previous cpDNA data [15], there was no sharing 
of haplotypes between China and Japan for either ITS 
(Figure la) or Tpi (Figure 2a). 

For each of these nuclear markers, the tree topologies 
recovered from MP and ML analyses were similar, and 
only the MP strict consensus trees are shown. According 
to the ITS tree (Figure 3), rooted with three outgroup 
species, P, arguta was monophyletic (with MP/ML 
bootstrap values of 100/100%) and consisted of two major 
clades corresponding to var. sinensis from East China 
(100/99%) and var. arguta from South Japan (84/69%). In 
addition to these lineages previously identified by cpDNA, 
the ITS marker revealed two novel subclades, 1 (98/89%) 
vs. 2 (99/98%), predominant in the southern vs. northern 
parts of South Japan (Kyushu, Shikoku vs. central Honshu) 
but with apparent range overlap in the Kii Peninsula of 
south-central Honshu (population J5; Figure la). As 
can be seen from the unrooted ITS haplotype network 
(Figure lb), the Chinese and Japanese haplotypes were 
separated from each other by 28 steps, and the two 
Japanese subclades by 13 steps. 

The Tpi phylogeny, rooted with Hydrangea chinensis 
(Figure 4), had lower resolution than the ITS tree but 
still recovered all Chinese haplotypes as monophyletic 
(85/87%), and the same applied to the majority of those 
from Japan (100/100%). Both clades, however, formed a 
trichotomy with a deviant clade (86/93%) of the two 
remaining haplotypes (T5, T15) from Japan. These latter 
haplotypes, which showed no distinct geographic distri- 
bution (Figure 2a), occupied an intermediate position in 
the unrooted Tpi haplotype network (Figure 2b), yet 
with somewhat closer relationships to China rather than 
Japan (four versus nine mutational steps apart from each 
group). 

For both ITS and Tpiy significant phylogeographic struc- 
ture was obvious at the species-level and within each of 
the two varieties (A/^st > G^st> all P values < 0.01; Additional 
file 6: Table S6). Hierarchical AM OVA revealed that most 
of the total nuclear sequence variation was partitioned be- 
tween them (ITS: 70.3%; Tpi: 66.7%; Table 1). Neverthe- 
less, for each variety {sinensis I arguta), Nst values were 
high at each nuclear sequence marker (ITS: 0.751/0.907; 
Tpi: 0.458/0.775) (Additional file 6: Table S6), reflecting 
that the majority of haplotypes were population specific 
(ITS: 32/33; Tpi: 26/31; see also Figures la and 2a). 

Nuclear microsatellites 

In 272 individuals from 14 populations of P. arguta, we 
detected a total of 220 alleles across 7 nSSR loci, with 17 
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100° 110° 120° 130° 140° 

Figure 1 Distribution of ITS ribotypes and the 95% plausible network of ribotypes in Platycrater arguta. (a) Distribution of ITS ribotypes. 

All ribotypes, except for H24 denoted by orange colour, are population-specific, and represented with different colours corresponding to the 

major phylogroups. Population codes are identified in Additional file 1: Table SI. (b) TCS-derived network of genealogical relationships between 

the 33 ribotypes. The small, solid black circles represent missing ribotypes. The sizes of circles are approximately proportional to sample size (n), 

with the smallest circles representing n = 1 and the largest representing n = 10. 
J 




100° 110° 120° 130° 140° 

Figure 2 Distribution of Tpi haplotypes and the 95% plausible network of Tpi haplotypes in Platycrater arguta. (a) Distribution of Tpi 
haplotypes. Population codes are identified in Additional file 1: Table SI. The distributions of shared haplotypes are denoted by colour, while 
private haplotypes are white, (b) TCS-derived network of genealogical relationships between the 31 haplotypes. The small, solid black circles 
represent missing haplotypes. The sizes of circles are approximately proportional to sample size (n), with the smallest circles representing n = 1 
and the largest representing n = 21. 
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Figure 3 Phylogenetic relationships of ITS ribotypes (H1-H33) of Platycrater arguta. Individuals of Hydmngeo anomola, H. chinensis and 
Schizophrognno hydrangeoides were used as outgroup taxa. Phylogenetic analyses using maximum parsimony (MP) and maximum likelihood (ML) 
produced trees with the same topology regarding major lineages. Only the MP strict consensus tree is presented. Numbers above and below the 
branches indicate, respectively, MP and ML bootstrap values (> 50%). 



to 48 alleles per locus. Average allele number (A/a) was 
higher in China (mean: 66) than in Japan (mean: 43), but 
allelic richness {Rs) and expected gene diversity {Hs) 
were very similar (China: Rs = 4.69, Hs = 0.734; Japan: 
Rs = 4.68, Hs = 0740; Additional file 1: Table SI). Within- 
population Fis values ranged from -0.054 to 0.480, with 
an average of 0.246. 

The STRUCTURE analyses provided strongest support 
for K = 7, both when considering the probability of the 
data LnP(D) and AK (Additional file 7: Figure SI). With 
K = 2, individuals from one population in China (C7) 
and Japanese populations formed a joint cluster, while 
the two varieties remained separate from K=3 upwards 
(Figure 5). At K= 7, five clusters were recovered in China 
(var. sinensis) and two in Japan (var. arguta), where most 
populations (J2-J5) comprised individuals representing 
both clusters (Figure 5). 

The overall Rst was 0.340, reflecting strong genetic 
differentiation over all populations and within each var- 
iety (var. sinensis I arguta: 0.273/0.366). However, in 
contrast to ITS and Tpi, the total genetic variance par- 
titioned between the two varieties was relatively low 
(13.5%; Table 1), suggesting pronounced allele sharing. 
Wilcoxon s test and the complementary mode-shift test 
provided little evidence for recent population bottle- 
necks, except for three populations from Japan (J2, J4, 



J6) and partly depending on the test applied (Additional 
file 8: Table S7). 



Isolation-with-migration (IM) analyses 

For this analysis, largest nonrecombining blocks of nDNA 
(ITS, Tpi) sequences were employed together with the 
nSSR markers and the previously obtained cpDNA se- 
quences [15]. The maximum-likelihood estimates (MLEs) 
and the 90% highest probability density (HPD) intervals of 
the six IMA-derived parameters are shown in Table 2, and 
their marginal posterior probability (MPP) distributions 
are illustrated in Additional file 9: Figure S2. Based on the 
geometric average mutation rate calculated (5.79 x 10"^ 
mutations per locus per year), these parameter estimates 
were converted to absolute values of years or individuals 
(Table 2). The split between var. sinensis and var. arguta 
was dated to about 889,358 yr BP, with a 90% HPD inter- 
val ranging from 509,438 to 1,193,295 yrBP (Additional 
file 9: Figure S2a, Table 2). For the ancestral (©a) and des- 
cendant (0c arid 0j) population rate parameters, both var. 
sinensis (Afc = 1.13 x 10^) and var. arguta (A/j = 6.00 x 10^) 
experienced a marked increase in effective population size 
(A/g) relative to that of their common ancestor (A/a = 
2.73 X 10^) (Additional file 9: Figure S2b, Table 2). Peak 
posterior estimates of post-divergence migration were 
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Figure 4 Phylogenetic relationships of Tpi haplotypes (T1-T31) 
of Platycrater arguta. Hydrangea chinensis was used as outgroup 
taxon. Phylogenetic analyses using maximum parsimony (MP) and 
maximum likelihood (ML) produced trees with the same topology 
regarding major lineages. Only the MP strict consensus tree is 
presented. Numbers above and below the branches indicate, 
respectively, MP and ML bootstrap values (> 50%). 



effectively zero in both directions (mc-j = 
Additional file 9: Figure S2c, Table 2). 



mj_c = 0.005; 



Bayesian skyline plot analysis of historical demography 

The EBSPs show that both varieties maintained low but 
stable effective population sizes (A/g) up to approximately 
0.4-0.5 Ma, and appear to have experienced a constant 
increase since (Additional file 10: Figure S3). While 
depicting a demographic trend over time, this analysis, 
however, cannot precisely estimate A/g because of the 
very broad confidence intervals (see Additional file 10: 
Figure S3). 



During the LIG (Figure 6b), the species' potential range 
in China was apparently reduced compared to the present. 
Also in the LGM (Figure 6c) only small and disjunct areas 
are predicted as suitable in the Wuyi/Yandang Mts. In 
contrast, in Japan, the species' current range is more or 
less similar to that during the LIG (Figure 6b), while dur- 
ing the LGM suitable habitat apparently contracted to the 
south (Kyushu, Shikoku) and the more northerly located 
Kii Peninsula (Figure 6c). Most strikingly, our LGM distri- 
bution model indicates almost no hospitable areas on the 
exposed ECS continental shelf, excepting those in the far 
east (i.e., extending from Kyushu to the delta region of the 
palaeo-Yellow River). 

Discussion 

Intraspecific divergence and post-divergence gene flow 

Our results provide strong evidence that the two recog- 
nized varieties of the WTD understorey shrub P. arguta in 
East China (var. sinensis) and South Japan (var. arguta) 
comprise reciprocally exclusive haplotypes for both ITS 
and Tpi (Figures 1 and 2). This is concordant with patterns 
for cpDNA [15] and consistent with the hypothesis that 
Chinese and Japanese conspecifics diverged over multiple 
glacial periods without inter-population gene flow. How- 
ever, relationships inferred from the STRUCTURE analysis 
of nSSRs with = 2 seem to conflict with this interpret- 
ation, as one Chinese population (C7) clusters with those 
from Japan (Figure 5). Given the congruence between nu- 
clear and chloroplast sequence data, the sharing of nSSR 
alleles most likely reflects homoplasy and/or incomplete 
lineage sorting rather than admixture through rare long- 
distance dispersal and/or migration across the ECS land 
bridge. In fact, our IMA analysis of combined multi-locus 
data (cpDNA, ITS, Tpi, nSSRs) recovered only close to 
zero signals of bidirectional post-divergence gene flow 
between the two varieties following their estimated di- 
vergence in the mid-Pleistocene, c. 0.89 Ma (90% HPD: 
1.2-0.51 Ma). Although this timing perfectly matches 
the respective point estimate drawn from cpDNA alone 
[15], single-locus cpDNA analysis suffered from a lack 
of convergence of the lineage divergence parameter {t), 
a problem that has been overcome with the present 
multi-locus approach (Additional file 9: Figure S2). 



Present and past ecological niche models 

The AUC value for the current potential distribution of 
P, arguta was 0.981, indicating a good predictive model 
performance. For the present (Figure 6a), these predicted 
areas mainly include the species' known distribution 
ranges in East China (Wuyi/Yandang Mts.) and South 
Japan (Kyushu, Shikoku, Kii Peninsula/south-central 
Honshu, and the Pacific Ocean side of central Honshu). 
Further suitable habitat is predicted in central-eastern 
Taiwan, but where the species is not known to occur. 



Broad-scale biogeographic history of R arguta 

According to Kimura [86,87], there were three main 
stages of land connections between the Eurasian mainland 
and the Japanese/Ryukyu Islands since the Late Miocene, 
with most of the ECS sea-floor exposed at about 7.0- 
5.0 Ma, 2.0-1.3 Ma, and 0.2-0.015 Ma, whereby the latter 
interval includes land bridge formations in both the penul- 
timate (Riss) and last (Wiirm) glacial periods [88-90]. Our 
estimated divergence time and its HPD intervals fall into 
the period (1.3) 1.0-0.2 Ma, the so-called 'Ryukyu Coral 
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Table 1 Hierarchical analysis of molecular variance (AMOVA) of ITS and Tpi sequences and nSSRs from 14 populations 
of Platycrater arguta var. sinensis (East China) and var. arguta (South Japan) 



Data sets/regional grouping Source of variation 
of populations 


d.f. 


Sum of squares 


Variance 
components 


Variance 
explained (%) 


F- or /?-stati sties (P) 


ITS 












var. sinensis (China) vs. var. arguta (Japan) 












Among groups 


1 


528.415 


13.954 


70.28 


FcT = 0.703** 


Among populations within groups 


12 


312.118 


4.912 


24.74 


Fsc = 0.833** 


Within populations 


58 


57.315 


0.988 


4.98 


FsT = 0.950** 


Tpi 












var. sinensis (China) vs. var. arguta (Japan) 












Among groups 


1 


357.190 


4.053 


66.70 


FcT = 0.667** 


Among populations within groups 


12 


175.449 


0.979 


16.11 


Fsc = 0.484** 


Within populations 


187 


195.305 


1.044 


17.19 


F5T = 0.828** 


nSSRs 












var. sinensis (China) vs. var. arguta (Japan) 












Among groups 


1 


25515.13 


96.12 


13.52 


/?CT = 0.135** 


Among populations within groups 


12 


7106.99 


1 78.93 


25.17 


Rsc = 0.290** 


Within populations 


530 


231038.92 


435.92 


61.31 


/?ST = 0.390** 



d.f., degrees of freedom; grouping of populations: Chinese group (populations C1-C7), Japanese group (populations J1-J7; for population codes see Additional file 
1 liable SI, Supporting information); P< 0.001. 

Estimators for Tpi and ITS were calculated based on the infinite alleles model (F-statistics) and those for nSSRs on the stepwise mutation model (/?-statistics). 



Sea Stage; when sea level rose tremendously and most 
formerly exposed land area submerged [87]. Based on 
these palaeo-data, and with no genetic indication of long- 
distance dispersal, the disjunct distribution of P. arguta 
across the ECS is implied to have resulted via vicariance, 
that is, the species was probably more widely distributed 



on the ECS land bridge in the early Pleistocene, around 
2.0-1.3 Ma, while divergence was subsequently driven 
by population contraction and extinction through land- 
bridge submergence. 

Unfortunately, it is not possible to test this 'early- 
Pleistocene expansion' hypothesis further using ENM, 




C1 C2 C3 C4 C5 C6 C7 J1 J2J3J4J5J6J7 



Platycrater y ax. sinensis (China) P. var. arguta (Japan) 

Figure 5 The proportion of genetic clusters detected by STRUCTURE analysis for the model with peaks at K= 2 and K= 7. The smallest 
vertical bar represents one individual. The assignment proportion of each individual into population clusters is shown along the y-axis. Note that 
STRUCTURE provided strongest support for K= 7, both when considering the probability of the data LnP(D) and AK (see text and Additional file 7: 
Figure SI). 
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Table 2 Maximum-likelihood estimates (MLEs) and lower and upper bounds of the 90% highest posterior density 
intervals (HPD90lo and HPDQOh./ respectively) of demographic parameters of Platycrater arguta from the IMA analysis 
of multi-locus data (cpDNA, ITS, Tpi, nSSRs) 



Estimates 


0c 


0j 


0A 


mc-j 


mj.c 


t 


Nc 


A/j 


A/a 


2NcMc-j 


2NjMj.c 


T (years BP) 


MLE 


106.596 


56.652 


25.730 


0.005 


0.005 


5.150 


113047 


60080 


27287 


0.00005 


0.00003 


889,358 


HPD90lo 


55.136 


29.992 


1 1 .027 


0.005 


0.005 


2.950 


58472 


31807 


11694 


0.00003 


0.00002 


509,438 


HPD90Hi 


231.570 


109.971 


69.839 


0.085 


0.085 


6.910 


245584 


116626 


74065 


0.00010 


0.00005 


1193,295 



Population rate parameters Go ©j/ and 0a refer to tine scaled effective population sizes (A/g) of var. sinensis (East China), var. arguta (South Japan), and the 
ancestral population, respectively, mc-j and mj.c are the scaled migration rates forward in time from var. sinensis to var. arguta and vice versa. Mc-j and Mj.c are 
the probabilities of migration from var. sinensis to var. arguta, per gene copy per generation and vice versa. 2NcMc-j and 2NjMj.c are the effective migration rates 
(number of migrants per generation), f is the time since ancestral population splitting in mutational units. 

All estimates include the per gene mutation rate u, which is equal to the geometric mean of the mutation rates of all the loci. 0c, 0j, 0a/ ^c-j/ '^j-o and f are 
scaled by the mutation rate, while Nq, A/j, A/a, 2NcMc-j, 2NjMj.c and T are scaled by individuals or years. 




Figure 6 Modelled climatically suitable areas of Platycrater arguta in East Asia at different times, (a) the present; (b) the Last Interglacial 
(LIG/Eemian: c. 130,000-114,000 yrBP); and (c) the Last Glacial Maximum (LGM: c. 21,000-18,000 yrBP). The current ecological niche model was 
established with six bioclimatic data layers on the basis of 59 sites of presence records of the species (black dots) using MAXENT 3.2.1 [76] and 
then projected onto a set of climatic variables simulated by MIROC 3.2 [80] to infer the extent of suitable habitats during the LGM and the LIG 
(see text). The map in (c) reflects changes in coastline and shelf exposition during the LGM due to lowered sea level (-110 m than at present; 
e.g., [7]). Yellow River and Yangtze Rive plaeo-channels in the exposed East China Sea basin are modified after [86]. The logistic value of habitat 
suitability is shown according to the grey-scale bars. 
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as no proxy- climate records are currently available for 
time periods earlier than the Last Interglacial (LIG) [91]. 
However, loess deposit and marine 6^^0 records indicate 
that both aridification and the influence of cold and dry 
winter monsoons in the northern subtropical areas of East 
Asia were less extensive before the mid-Pleistocene cli- 
mate transition (MTP), a 0.85 Ma [92]. Afterwards, glo- 
bal ice volume drastically increased along with a change 
in the dominant orbital cycle from 41,000 to 100,000 years, 
resulting in drier and colder climates during glacial pe- 
riods [93]. Given that P. arguta is a moisture-dependent 
and drought-sensitive shrub, early-Pleistocene environ- 
mental conditions may have still favoured a more con- 
tinuous distribution on the ECS land bridge, whereas 
subsequent glacials would have prevented the migration 
of individuals between China and Japan. Indeed, our 
LGM distribution model for P, arguta essentially shows 
no suitable habitat areas on the exposed ECS basin, with 
the exception of the delta region of the palaeo-YeUow 
River (Figure 6c). We therefore conclude that the ECS 
land bridge acted as a formidable barrier to the dispersal 
of P, arguta during the LGM and earlier cold periods of 
the Late Pleistocene. However, additional factors, other 
than climate-related niche requirements, may have had 
a role in preventing, or at least hampering such disper- 
sal. This is because similar (i.e., genetic and ENM) evi- 
dence suggests that the ECS land bridge served as an 
intermittent route of range expansion for most of the 
Late Pleistocene to counteract isolation between Chin- 
ese and Japanese populations of two widespread WTD 
forest tree species, namely Cercidiphyllum japonicum 
[12] and Kalopanax septemlobus [13]. Both are tall canopy 
trees with long generation times, large effective population 
sizes and high seed production [94,95], and additional 
traits generally considered important for population 
recovery, such as a generalist (i.e., wind) pollination/ 
dispersal system (C. japonicum) and vegetative repro- 
duction (C. japonicum, K, septemlobus). Similar traits are 
largely absent or only moderately developed in P, arguta 
(except for its wind-dispersed seed); the same applies to 
other rare and narrowly distributed understorey plants of 
the same forest biome with likewise ancient (early-to-mid 
Pleistocene) genetic breaks across the ECS (Croomia 
japonica: [14]; Kirengeshoma palmata: [16]; Ligularia 
hodgsonii: [17]). Moreover, in forest habitats particu- 
larly, one should not dismiss the effects related to 
edges (e.g., abundance to animal pollinators and/or dis- 
persers), which may either increase or decrease seed pro- 
duction and recruitment [96-98]. Hence, species in the 
forest interior, such as P. arguta and other understorey 
plants, should be negatively affected by fragmentation; 
by contrast, species using edge or transitional habitats, 
such as the riparian-dwelling C. japonicum and the 
semi-invasive K. semptemlobus, may be favoured by 



fragmentation [96,97]. We therefore propose that (i) 
range expansion in response to the formation of the gla- 
cially exposed ECS land bridge is species-specific; and 
(ii) predictions on the effects of this ephemeral land 
bridge, as corridor' or 'filter', have to account not only 
for habitat preferences per se but also for other bio- 
logical features of each species, especially its recruit- 
ment properties. 

Inferences of historical demography and range dynamics 
in P. arguta 

Our results provide strong evidence that both var. sinensis 
and var. arguta underwent past population growth follow- 
ing their inferred sundering in the mid-Pleistocene. This is 
supported by significantly positive R2 statistics for both 
ITS and Tpi (Additional file 5: Table S5), and is also evi- 
dent from the joint analysis of the four genetic data sets 
(cpDNA, ITS, Tpi, nSSRs) with IMA, suggesting a some- 
what larger effective population size (A/g) in each variety 
compared to their ancestral population (Additional file 9: 
Figure S2b, Table 2). We caution, however, that this esti- 
mate of ancestral A/g probably reflects post-vicariant 
conditions, and thus is biased downward, as complex 
population dynamics within the ancestral population, 
such as contractions and extinctions, are not accounted 
for by IMA [99]. Nonetheless, analyses of the combined 
nuclear and plastid sequence data with EBSPs were consist- 
ent in showing a strong increase of A/g in each variety from 
about 0.4-0.5 Ma onwards (Additional file 10: Figure S3). 
This is very similar to the point estimates of expansion 
times inferred from mismatch analyses of cpDNA alone, 
that is, c. 0.43 and 0.45 Ma for var. sinensis and var. 
arguta, respectively [15]. Hence, the present results sup- 
port our previous notion that this near-synchronized 
population growth may have been triggered by climate 
change at the beginning of Chinas 'Penultimate Intergla- 
cial Period' (c. 0.46-0.33 Ma; [100]). This also accords 
with palaeo-climate data, suggesting that during intergla- 
cials since the mid-late Pleistocene (c. 1.0-0.78 Ma) the 
warm, wet East Asian summer monsoons have intensified 
[101-103]. In consequence, this may have also led to an in- 
crease of suitable habitats for P, arguta throughout the 
warmer periods of the Late Quaternary. Partly consist- 
ent with this, the ENM shows larger and more contigu- 
ous distributions at least in northern South Japan [i.e. 
Kii Peninsula and (south-)central Honshu] at both the 
LIG and the present compared to the LGM (Figure 6). 
Although suitable habitats may have shrunk in East 
China at the LIG (Figure 6b), and in both regions during 
the LGM (Figure 6c), we found no direct (EBSP) evidence 
for population decline (Additional file 10: Figure S3) and 
recent population bottlenecks (Additional file 1: Tables SI; 
Additional file 5: Tables S5). Nonetheless, the ENM sug- 
gests that potential glacial refugia in East China were more 
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strongly affected by small-scale fragmentation compared 
to the situation in South Japan (Figure 6c) in general, 
and its southern parts (Shikoku, Kyushu) in particular, 
where such refugia might have even existed in nearby 
shelf areas. This inter-regional difference in glacial habi- 
tat suitability most likely reflects differences in topog- 
raphy and (micro-) climate, and should have differing 
consequences for population genetic structure and evo- 
lutionary history (see below). 

Contrasting Late Quaternary evolutionary histories of 
P. arguta in China and Japan 

In var. sinensis, there is a marked subdivision between 
populations in nuclear genes (ITS, Tpi, nSSRs) (Table 1), 
along with high haplotypic and allelic diversity (Additional 
file 1: Table SI), and the same holds true for cpDNA [15]. 
Together with our STRUCTURE result of five principal 
nSSR gene pools in this variety (Figure 5), all of these data 
indicate long-term population persistence and isolation 
over multiple glacial/ interglacial cycles. With no indication 
of latitudinal range shifts, past population growth in var. 
sinensis (see above) may also reflect, at least in part, re- 
peated downhill shifts in elevation range during glacials, 
perhaps by tracking favourable humidity conditions as im- 
posed by the East Asian monsoon in areas of high relief 
[104]. Over climatic cycling, such contiguous but localized 
range shifts along elevation gradients would have mini- 
mized bottlenecking and loss of genetic diversity, ultim- 
ately resulting in strong population subdivision [105]. 

In Japan, patterns of cpDNA diversity in var. arguta 
were previously interpreted to indicate southward retreat 
during glacials (to Kyushu and Shikoku) followed by 
expansion northward (up to central Honshu) during 
inter-Zpost-glacials [15]. However, the current nuclear 
diversity and ENM data cast doubt on the validity of 
the 'expansion-contraction' (EC) model for var. arguta. 
For example, nuclear genetic diversity within var. arguta 
populations is more or less evenly spread throughout 
the distribution (Additional file 1: Table SI), and neither 
haplotypic (/zs) diversity (ITS, Tpi) or allelic richness 
{R^; nSSRs) show a significantly negative relationship 
with latitude {P = 0.33-0.58; data not shown), as would 
be expected under a scenario of south-to-north (re-) 
colonization [106,107]. In addition, the ENM for the 
LGM indicates suitable habitats not only in the south 
(Kyushu, Shikoku, and nearby shelf areas) but also, 
more disjunctively, in the Kii Peninsula of south-central 
Honshu (Figure 6c). By contrast, large areas of contiguous 
suitable habitats were likely present at the LIG (Figure 6b), 
suggesting that opportunities for admixture occurred re- 
peatedly during inter-Zpostglacial times, unless hampered 
by the spread of evergreen forest [1,4]. 

Both of these latter predictions are basically bone out 
by the current nuclear data, even though each marker 



suggests a somewhat different historical scenario. Thus, 
the identification of two separate ITS lineages in Kyushu/ 
Shikoku vs. central Honshu, with an apparent overlap of 
range in the Kii Peninsula (pop. J5; Figures la and 3), is 
strongly suggestive of a two-refugia scenario and a narrow 
contact zone between the two lineages. Interestingly, this 
putative contact region coincides with a well-known inter- 
and intraspecific suture zone of Japans temperate fiora 
and fauna [14,108,109,12]. On the other hand, the Tpi data 
are more consistent with a multiple-refugia scenario, that 
is, haplotypes are largely restricted to particular (sub)re- 
gions and/or populations, while still showing a relatively 
small amount of shared polymorphisms, mainly (but not 
exclusively) among adjacent sites (Figure 2a). Finally, the 
nSSRs reveal extensive admixture in most populations 
(J2-5), excepting those in Kyushu (Jl) and central Honshu 
(J6, J7; Figure 5), suggesting that, in marked contrast to 
ITS, the contact area is much wider towards the south. 

This discordance among patterns of genetic structure 
observed in the four genetic data sets makes it difficult 
to specify exactly which historical processes occurred in 
var. arguta. This discordance (even among nuclear genes 
and relative to cpDNA) is not readily explicable by a single 
difference in marker properties [e.g., mutation rates, trans- 
mission genetics, effective population sizes, concerted evo- 
lution (only ITS)], but likely results from a combination of 
marker features, with both incomplete lineage sorting and 
secondary admixture acting at different temporal and 
spatial scales [83]. Taken together, we suggest that the gen- 
etic and ENM patterns found in var. arguta refiect a rela- 
tively ancient north-south vicariant event during glacials 
(sundering populations in central Honshu), followed by 
more recent climate -induced cycles of range contraction 
and expansion/admixture, with the latter producing more 
shallow patterns of divergence in more southerly areas 
(Kyushu, Shikoku, Kii). 

Conclusions 

This is the first study providing firm evidence that the 
ECS land bridge acted as a 'filter' during the last glacials 
in selectively preventing the dispersal of certain plant 
species of WTD forest, such as P, arguta, even though a 
near continuous belt of this forest biome is thought to 
have covered this land bridge during the LGM [2]. Our 
data emphasize the species-specific recruitment and 
range expansion response of WTD forest- dwellers to the 
formation of the ECS land bridge, and caution against 
an uncritical interpretation of reconstructed palaeo-forest 
biome data as guides of past range dynamics of constitu- 
ent species populations [8]. 

In addition, this multi-locus study strongly supports 
the two varieties of P, arguta in East China (var. sinensis) 
and South Japan (var. arguta) as genetically distinct units 
that evolved in strict allopatry since the mid-Pleistocene, 
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which concurs with previous findings inferred from 
cpDNA alone [15]. However, the combination of genetic 
structures from both nuclear and cytoplasmic loci can 
better depict the history of populations, demonstrating 
that (i) each lineage has undergone refugial isolation 
and divergence; and (ii) var. arguta likely experienced 
post-glacial admixture across a well-known suture zone. 
Yet, the extent of the species' overall distribution does 
not seem to have greatly changed over the last glacial- 
interglacial cycles. 
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effective population sizes (0) of var. sinensis (0c)/ var. arguta (0j), ar^d the 
ancestral population (0a). (c) The scaled migration rates forward in time 
from East China to South Japan (mcJ, and vice versa (/t]j_c). 

Additional file 10: Figure S3. Extended Bayesian Skyline Plots (EBSPs), 
inferred from cpDNA and nDNA (ITS, Tpi) sequence variation and 
depicting changes in effective population size (A/e) as a function of time 
for (a) Platycrater arguta var. sinensis (East China) and (b) var. arguta 
(South Japan). The thick solid black line is the median estimate, and the 
area delimited by the upper and lower grey lines represents the HPD 
95% confidence intervals for N^. 
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